#include "Dwarf.Unmanaged.Serial.h"

//Init arrays.
Solver :: Solver()
{
    lgrid = 0;
    rgrid = 0;
    result = 0;

    msize = 0;
}

//Dispose arrays
Solver :: ~Solver()
{    
    for (int i = 0; i < msize; i ++) 
    {
        delete lgrid[i];
        delete rgrid[i];
        delete result[i];
    }

    delete lgrid;
    delete rgrid;
    delete result;
}

static void Transpose(double** matrix, double** transpose, int n)
{
    for (int i = 0; i < n; i++)
    {
        transpose[i][i] = matrix[i][i];
        for (int j = i + 1; j < n; j++)
        {
            transpose[i][j] = matrix[j][i];
            transpose[j][i] = matrix[i][j];
        }
    }
}

// Serial based method for matrix-matrix multiplication.
void Solver::solve() 
{ 
    // transpose the right hand array to improve memory accessing
    double** rtransposed = new double*[msize];
    for (int i = 0; i < msize; i++)
    {
        rtransposed[i] = new double[msize];
    }
	Transpose(rgrid, rtransposed, msize);

    //Multiple lgrid by rgrid (transposed)
    //    Note that row and cell are factored out
    //    since the optimizers shouldn't hoist all of 
    //    the common addressing subexpressions without doing it. 
    //    (Strickly speaking it would be incorrect to do.)

    for (int i = 0; i < msize; i ++ )                                   // Loop for number of rows in left matrix.
    {
        double* row = result[i];
        for (int j = 0; j < msize; j ++ )                               // Loop for number of columns in right matrix.
        {
            double cell = 0.0;
            double* leftRow = lgrid[i];
            double* rightColumn = rtransposed[j];
            for (int k = 0; k < msize; k ++ )                           // Loop for number of columns in left matrix and rows in right matrix.
            {
				cell += leftRow[k] * rightColumn[k];

                // To see the effect of not transposing 
                // use the following as the inner loop
                //cell += leftRow[k] * rgrid[k][j];
                // To see the effect of the compiler/runtime optimizations,
                // comment out cell, row, and leftRow above and below.
                // Use the following as the inner loop
                //result[i][j] = result[i][j] + lgrid[i][k] * rgrid[k][j];
            }
			row[j] = cell;
        }
    }

	for(int i=0; i<msize; i++)
	{
		delete rtransposed[i];
	}
	delete rtransposed;
}